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We study the interplay of population growth and evolutionary dynamics using a stochastic model 
based on birth and death events. In contrast to the common assumption of an independent pop- 
ulation size, evolution can be strongly affected by population dynamics in general. Especially for 
fast reproducing microbes which are subject to selection, both types of dynamics are often closely 
intertwined. We illustrate this by considering different growth scenarios. Depending on whether 
microbes die or stop to reproduce (dormancy), qualitatively different behaviors emerge. For co- 
operating bacteria, a permanent increase of costly cooperation can occur. Even if not permanent, 
cooperation can still increase transiently due to demographic fluctuations. We validate our analysis 
via stochastic simulations and analytic calculations. In particular, we derive a condition for an 
increase in the level of cooperation. 
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I. INTRODUCTION 

The time evolution of size and internal composition 
of a population are both driven by discrete birth and 
death events. As a consequence, population dynamics 
and internal evolutionary dynamics are intricately linked. 
The biological significance of this coupling has previ- 
ously been emphasized [THS]- Those studies mostly em- 
ploy density-dependent fitness functions to phenomeno- 
logically derive sets of coupled deterministic equations for 
the size and composition of populations in various eco- 
logical contexts. While those studies correctly describe 
the evolutionary dynamics of large population sizes, they 
do not account for stochastic effects arising at low pop- 
ulation sizes. These demographic fluctuations are natu- 
rally described in the theoretical framework of stochastic 
processes based on elementary birth and death events as 
recently introduced [10]. In particular, this approach al- 
lows one to explore the role of fluctuations in populations 
with a time-varying population size. 

To understand such interdependence of population and 
evolutionary dynamics, it is instructive to first review the 
decoupled and deterministic formulations of both. Evo- 
lutionary game theory is a well defined framework to de- 
scribe the temporal development of different interacting 
traits or strategies [11] [12] . It has been established as 
a standard approach to describe evolutionary dynamics 
if the fitness is frequency-dependent, i.e. if the fitness of 
a certain strategy depends on the abundance of other 
strategies within the population. Within the most ba- 
sic setup, well-mixed populations are assumed and the 
evolution of strategies is solely determined by fitness ad- 
vantages. The temporal development of the abundance 
Xs oi a. trait S follows a replicator dynamics |llfll3| . 

dtXs = (0s - 0) Xs. (1) 

A trait's abundance increases if its fitness (ps exceeds 
the average fitness 4> of the population. The frequency 



dependence, with (ps a function of the abundances x of 
all strategies, provoke non-linearities in Eq. ([T]). Start- 
ing from this standard approach, many specific exam- 
ples and extensions thereof have been studied |f21ll4j . 
This comprises, for example, the prisoner's dilemma, the 
snowdrift-game and other games in well-mixed popula- 
tions [11 13, 15 . ft further ranges from the role of spatial 
arrangements and network interactions |16H25j via cyclic 
dominance [71 l26ti33] . structured populations [34l l35] . 
modified update-rules [5S1 |3Z], multi-player games [35] 
and evolutionary algorithms [39' to the influence of inter- 
nal and external fluctuations [40H45j . While these mod- 
els consider a wide range of evolutionary aspects, they 
mostly rely on one key assumption, a decoupled, con- 
stant population size. 

In contrast, population dynamics focuses on the time 
evolution of the population size and how it is determined 
by environmental impacts like limited resources or sea- 
sonal variations. The dynamics is typically described by 
differential equations of the form [ISHIS] 

dtN^F{N-t), (2) 

where T{N\t) may explicitly depend on time [46]. The 
most prominent example is logistic growth [33]. While a 
small population grows exponentially, the growth rate de- 
creases with increasing population size due to limitations 
of resources and the population size is bounded below a 
maximum carrying capacity. 

Illustrative examples of dynamical changes in the 
population size comprise bacterial and other microbial 
populations [5Dl - f5^ : A surplus in nutrients or other 
metabolism related factors, can lead to an immediate 
and strong growth of the population while resource lim- 
itations or antibiotics and other detrimental factors can 
imply a stop in growth or even an abrupt death of single 
individuals. Even for only slightly varying environmen- 
tal conditions, a fixed population size is thus rather the 
exception than the rule. 
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But microbes not only show rich population dynamics, 
they are also subject to diverse evolutionary forces [53l - 
[57] • Microbes live in interacting collectives of differ- 
ent traits. Evolution is ubiquitous and strong forms 
of frequency-dependence can be observed. Public good 
scenarios where a metabolically costly biochemical prod- 
uct is shared among individuals are of particular interest 
from an evolutionary perspective, see e.g. [Slj.SS, 58-60]. 
This includes, for example, nutrient uptake, like disac- 
charides in yeast [GTMBBj . collective fruiting body forma- 
tion [M] [5S] , or the active formation of biofilms [5ll [571 
I66| 167] . An example regarding iron uptake is considered 
below in more detail (68n70) . Furthermore, synthetical 
microbial systems have been considered [711 [72] ■ 

Motivated by these recent studies of microbial systems, 
we here investigate the consequences of such an interde- 
pendence between evolutionary and population dynam- 
ics. Employing a previously introduced theoretical ap- 
proach 10. , we study the influence of different growth 
scenarios in combination with demographic fluctuations. 

The outline of this article is the following. In Scction|TT] 
we discuss the stochastic dynamics and its deterministic 
approximation. Furthermore, we consider the limits in 
which the model maps to standard (deterministic and 
stochastic) formulations of evolutionary dynamics. In 



trait A 
trait B 



Section III we consider the dilemma of cooperation in 
growing populations. Here, an increase of cooperation 
can be observed which is analyzed in detail. In partic- 
ular, we discuss the outcomes for two different growth 
scenarios, i.e. a reproduction-dynamics which either is 
balanced by death events or simply arrests in the sta- 
tionary case. Finally, we close with a short conclusion in 
Section [TV] 
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FIG. 1: (Color online) The per capita birth and death rates 
for two different traits, A [light gray (red)] and B [dark gray 
(blue)]. Each rate depends on a global, trait-independent 
and a relative, trait-dependent part. While the global and 
relative fitness terms, g and /a/b affect the birth rates, the 
global and relative weakness terms d and wa/b determine the 
death rates. 



We write 

Gs^g{S-,N)fs{S), Ds = d{x,N)ws{S!), (4) 

and refer to g{x, N) and (i(x, N) as global birth-fitness 
and global weakness, respectively. The trait-dependent 
terms are the relative birth-fitness fs{x) and the relative 
weakness 'ws{x)^U^. While birth- fitness terms affect the 
birth rates, weakness terms determine the expected sur- 
vival times of individuals and hence their viability. A 
short illustration of the stochastic processes is given in 
Fig. |l] for the case of two different traits. 

To specify the relative fitness terms, we follow the stan- 
dard approach of evolutionary game theory [llj . and as- 
sume them to depend linearly on the frequencies x. Let 
V be the payoff matrix for birth events. Then, the cor- 
responding fitness vector for all traits is defined as 



= 1 + sVx. 



II. COUPLING OF EVOLUTIONARY AND 
POPULATION DYNAMICS 

A. Microscopic Model 

We consider a population of M different traits. Each 
trait S is represented by Ns individuals, such that the 
state of the population is given by TV = {Ni,N2, Nm)- 
We further denote the frequencies of all different traits by 
X ~ N/N with N — J2s being the total population 
size. The stochastic evolutionary dynamics is formulated 
in terms of per capita birth and death rates, Gs and Ds, 
respectively. The total rate for the abundance of trait S 
to increase or decrease by one individual is given by 



GsNs, 



Fc 



DsNs. 



(3) 



The various biological factors determining each rate can 
be split up into two parts, a global and a relative contri- 
bution. While the global term is trait-independent and 
affects all traits in the same manner the relative term is 
trait-dependent and sets the differences between traits. 



(5) 



Following standard formulations, the selection strength s 
defines the relative weight of a frequency-dependent part 
with respect to a background-fitness set to 1 ^01 [33]. 
As will become clear in the following, it is convenient to 
make use of normalized fitness values. 



(6) 



where (p — J2s 4>sXs such that / — fs^s — 1- With- 
out loss of generality, this choice separates global and 
relative parts in such a way that the dynamics of the pop- 
ulation size depends only on the global functions g and d; 
see also the following Eq. (8a). An analogous approach 



with a payoff matrix V for death events can be used to 
obtain the frequency-dependent weakness functions ws, 
which are also taken as normalized, "^gWsxs = 1- Of 
course, a more general, non-linear frequency-dependance 
for both relative functions can readily be taken into ac- 
count. For example, in microbial systems the fitness of 
an individual or the whole community depends in an 
intricate way on a plethora of factors, e.g. the abun- 
dance of individuals, secretion and detection of signaling 
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molecules, toxin secretion leading to inter-strain compe- 
tition and changes in environmental conditions. Non- 
linear frequency-dependent fitness-functions might help 
to account for such factors, see e.g. [551IS^ . 

In general, the global terms g{x, N) and d{x, N) 
depend on the population size and are frequency- 
dependent. Limited growth is one example of size- 
dependence. In such a setting, small populations start 
to grow exponentially but growth is bounded due to lim- 
ited resources, e.g. d{x, N) increases with N. Frequency- 
dependent terms can, for example, occur in public good 
situations, as discussed in Section IIIII 



B. Dynamics 

The per-capita birth and death rates, Eqs. Q, define a 
continuous-time Markov process [731 IZl] • It is described 
by a master equation for the probability density P{N;t) 
to find the population in state N at time t: 



dP{N;t) 
Jt 



= J2 - ^)GsNs 
s 

+ {E+ -l)DsNs]P{N;t). (7) 



Here, Eg are step operators increasing/decreasing the 
number of individuals of trait S by one [74], e.g. 



PiNi,...,Ns±l,...,NM;t). 



For a reference it is instructive to first consider a de- 
terministic limit where both fluctuations and correlations 
can be neglected. Then, upon factorizing higher mo- 
ments of the probability density [731 El], one finds a 
closed set of equations for the expected frequencies xs 
and the total population size TV: 

dtN^[g (f , N)f~d {x, N)w]N, (8a) 
dtxs = g{x, N) [fs{x) - / ] xs 

— d{x, N) [ws{x) — w\xs 1 (8b) 

where J — w — 1 according to Eq. ([6|. To unclutter 
notation, we have not explicitly marked the expectation 
values in Eqs. ([S]) but use the same notation as for the 
stochastic variables. 

This set of coupled non-linear equations resembles 
other deterministic approaches [iHll [SHS] and has a sim- 
ple interpretation. Eq. (8a) describes the population dy- 
namics. As is typical for a deterministic approach, the 
dynamics does not depend on the global birth-fitness, 
and the global weakness, d, separately, but only on their 
difference. Eq. (8b) describes the internal evolution of 



the population: The time evolution of the frequency of 
a strategy S is given by the interplay between a growth 
and a death term. Each of them consists of a relative 
term measuring the surplus of the fitness/weakness rela- 
tive to the corresponding population average. The weight 



of these terms are given by the respective global fitness 
functions, g and d. During phases of population growth, 
where g > d holds (see Eq.([8a[)), the growth term and 
hence differences in relative birth fitness dominate the 
internal evolution of the population. Similarly, weakness 
differences are the main evolutionary driving forces dur- 
ing population decline. 

From these considerations it follows that both the time 
scale of population and evolutionary dynamics have a 
crucial impact on the dynamics. This is obvious if the 
time-scales are similar. Such biological situations have 
been observed in many examples, see e.g. [75H75] . But 
also if evolution happens on longer time-scales than ecol- 
ogy this coupling can affect the evolutionary outcome as 
we show in the following. 

Importantly, fluctuation cannot be ignored in gen- 
eral but can change evolutionary dynamics dramatically. 
Then, the deterministic approach given by Eqs. ([8| is 
not adequate. This regards for example flxation and ex- 
tinction events but also the evolution of first and higher 
moments of a trait's abundance. For a proper descrip- 
tion, one has to take the full stochastic dynamics and 
master equation ([t]) into account. One example, where 
fiuctuations drastically change the outcome is given in 
the following Section [TTTl 



C. Mapping to Standard Approaches: Replicator 
Dynamics and the Moran Process 

We now consider in which limits and to what extent our 
stochastic approach resembles the standard approaches 
of evolutionary dynamics. Let us first consider the spe- 
cial case where the global rates g{x, N) = g{N) and 
d{x, N) = d(N) are frequency-independent and the ensu- 
ing deterministic dynamics exhibits a stable fixed point 
N* in the population size. Then, birth and death events 
exactly balance each other, g{N*) = d{N*), such that N* 
is fixed, dtN* = 0. This is, for example, the case if the 
population size evolves according to a logistic growth law 
and the carrying capacity has been reached. In the de- 
terministic limit, the internal dynamics, Eq. ([8])(b), sim- 
plifies to 

dtXs - g{N*) [fsix) - f-wsix) + w] xs. (9) 

The fraction xs evolves like in a standard replicator 
equation, similar to Eq.([I]). It is the difference of both 
relative terms, the effective fitness fs — ws, which deter- 
mines internal evolution. Compared to Eq. ([T]), the ad- 
ditional constant prefactor g{N*) in Eq. (|9| just rescales 
the time-scale on which internal evolution occurs |41j . 

Furthermore, also the full stochastic formulations of 
our model and the standard stochastic approaches with a 
fixed population size resemble each other. In those stan- 
dard approaches, the birth of one individual is directly 
coupled to the death of another one. The dynamics is 
described by update rules. For example, for the time- 
continuous formulation used here, the stochastic dynam- 
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ics can be described by the Moran process [101 STJ |331 1^ 
FfOHSTj |101j . In our formulation, this process holds in the 
limit where the fixed point of the population size, TV*, 
is linearly stable with a large stability coefficient |102j . 
Then, a birth event is directly followed by a death event 
and vice versa. The effective rate for such a combined 
birth-death event is given by, 

Ts-yS' = ^S'^2S'^S-f<6 + rs^0rs'-i.2S'- (10) 

The strength of fluctuations in the fraction of a certain 
species is of the order l/\/N* and the transition rate 
Vs^S' follows by the logic of an urn-model where, fltness- 
dependent, individuals reproduce to substitute other, 
randomly chosen, individuals [iP l |IT1 liH l [75 1 [50 ] . 

Beyond the Moran process, however, if N* is not lin- 
early stable with sufficiently high stability coefficients, 
then birth and death events do not strictly follow each 
other. Depending on the stability of the fixed point, evo- 
lutionary paths deviating from A'^* by more than one indi- 
vidual have to be taken into account to derive an effective 
rate for a combined birth-death event. 

In general, the population size changes with time, 
N = N{t). For frequency-independent global rates, the 
deterministic limit of the internal evolutionary dynamics 
resembles the form of a replicator equation, 

dtN = [g{N)-d{N)]N, (11a) 
dtXs={9{N)[fs{S)^J] (lib) 
-d{N)[ws{S)~w])xs. (11c) 

However, in contrast to Eq. ([T]), both relative fitness 
terms, / and w, are now weighted by the global rates. 
This has important implications. While in growth phases 
with g > d the relative birth fitness fs dominates the 
dynamics, the relative weakness functions ws dominate 
during population-decline, g < d. Moreover, the time- 
varying population size also leads to a changing strength 
of fluctuations ^ 1/ ^N{t). In particular, when fitness 
differences are weak and the dynamics is close to neu- 
tral evolution, such a change might have strong conse- 
quences mi HI ESUHH IH3] ■ 

III. THE DILEMMA OF COOPERATION IN 
GROWING POPULATIONS 



load of production, permanently run the risk to be under- 
mined by non-producing free-riding strains. An example 
is provided by the proteobacterium Pseudomonas aerug- 
inosa [BS'-TfT. To facilitate the metabolically important 
iron-uptake, these microbes produce siderophores which 
they release into the environment. Given the high bind- 
ing affinity to iron, these proteins are capable of scav- 
enging single iron atoms from larger iron clusters. The 
iron-siderophore complex can then be taken up by the 
bacteria, ensuring their iron supply. However, as every 
bacterium, not only the producing ones, can take ad- 
vantage of the released siderophores there is a dilemma 
of cooperation: While it would be optimal for the whole 
population to cooperate, cooperators are endangered due 
to their reproduction disadvantage. 

In addition to the evolutionary dynamics, microbial 
colonies are also subject to strong changes in popula- 
tion size [SHHS^ [M] . While in the presence of nutrients, 
small colonies grow exponentially, growth is bounded 
due to limitations in resources or deteriorating environ- 
mental conditions. This includes insufficient amounts of 
nutrients, a lack of oxygen or a poisoning by metabo- 
lites. Eventually the colony size remains constant or 
even declines again [50] . Given by the exact interplay 
of these detrimental and other environmental factors, 
and differing from species to species, growth dynamics 
varies between two scenarios }851 186] . First, bacteria can 
switch into a dormant state where individuals stay alive 
but regulate reproduction rates and metabolic activity 
towards zero {dormancy scenario). Depending on en- 
vironmental conditions dormancy can increase survival 
chances. For example, in the presence of antibiotics, this 
downgraded metabolism can make bacteria less vulner- 
able leading to persistence |57H5Uj . or dormancy might 
hedge a population against strongly fluctuation environ- 
ments [SGllQOlEI]. Second, environmental conditions can 
lead to death rates increasing with the population size N 
while birth rates are only slightly affected [22] . The popu- 
lation, therefore, reaches a state of dynamical maintained 
population size with the death rates exactly balancing the 
birth rates {scenario of balanced growth). In many popu- 
lations, a situation in-between both scenarios is observed. 
In pathogenes like P. aeruginosa, the fraction of individ- 
uals transferring to the dormancy state varies between 
20% and 80% [S^ . In the following we consider both sce- 
narios and their impact on internal evolution separately. 



To exemplify the importance of coupling and fluc- 
tuations offered by our approach, we here study the 
dilemma of cooperation in growing populations. This 
is motivated by the dynamics observed in microbial 
biofilms where strong forms of cooperation can be ob- 
served [5Tl|S51|571i59, 60, 67 . Single individuals produce 
metabolically costly products which they release into the 
environment to support, for example, biofilm formation 
or nutrient depletion. As these products are available for 
other bacteria in the colony, the cooperating individuals 
are producers of a public good, and, by having the extra 



A. The Balanced Growth Scenario 

Let us first study the balanced growth dynamics where, 
in the stationary state, birth and death events are both 
present, but exactly balance each other such that the 
population size is about constant. We consider a popu- 
lation which consists of two traits, cooperators (C) and 
free-riders {F). The total number of individuals in the 
population is given by iV = Nc + Np and the fraction 
of cooperators hy x = xq — Nc/N. The relative birth 
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fitness, fs (0s, if not normalized), accounts for the re- 
production disadvantage of cooperating individuals. We 
study the well-know prisoner's dilemma [TT] [103] : 



be —I + s{bx — c), 
tip =1 + sbx, 
(f> =1 -1- s{b — c)x. 



(12) 



As introduced in Section |TTj the frequency-dependent 
part is weighted with the strength of selection s. Indi- 
viduals obtain a benefit b from direct interaction with 
cooperators, while only cooperating individuals have to 
pay the cost c for producing the public good. For the 
resulting normalized fitness functions, fg = (f>s/4'i 
inequality fc < fp always holds; within the same pop- 
ulation, the reproduction rate of cooperators is always 
smaller than the one of free-riders. 

In the following, we take the payoff parameters to be 
constant, c — I and 5 = 3. Then, s directly sets the time 
scale of the internal evolution. The relative weakness 
is assumed to be trait- independent and constant, wc — 
wp = 1; free-riders and cooperators have equal survival 
chances. 

Furthermore, because cooperators are the producers of 
a public good, the overall growth condition of a popula- 
tion improves with a higher level of cooperation. We here 
choose the global birth fitness to increase linearly with 
the level of cooperation. 



g{x) = 14- px. 



(13) 



The parameter p scales the positive impact of the pres- 
ence of public good on the population. In the scenario 
of balanced growth, we consider death rates increasing 
with the population size. For specificity, we assume lo- 
gistic growth [49 and set 



d{N) = N/K. 



(14) 



K scales the maximal size a population can reach (car- 
rying capacity) as discussed in detail below. 

The master equation ([t]) describing the full stochastic 
dynamics then takes the form 



dPjNcNp) 
dt 



[iEc-l)9fcNc + iiEp-l)gfpNp 

+ {E+-l)dNc + (E'^-l)dNc] X 
P{Nc,Np). (15) 



To explore the dynamics, we performed extensive 
stochastic simulations. They were obtained by simulat- 
ing i = 1, . . . , i? different realizations with the Gillespie 



algorithm [§(4j, according to the master equation ( 15 1. In 
Fig. [2j we show the ensemble averages of the population 
size (A'') and the fraction of cooperators (x) given by 

(TV) =^A^,(i)/i?, (16a) 

i 

(x)=^^c.(i)/E^'W- (16b) 



This choice for the average naturally accounts for the 
fact that realizations with a larger populations size have a 
larger weight. It is especially important for biological sit- 
uations where several realizations exist at the same time, 
e.g. |72j . In such an ensemble cooperation can increase in 
principle if there is a positive correlation between popu- 
lation size and the fraction of cooperators. The existence 
of this effect, also known as Simpson's paradox, has been 
shown recently by Chuang et al. for microbial popula- 
tions |72j . Here we want to understand the dynamics 
underlying this correlation underlying cooperation. 

Starting with a small population, the system size grows 
exponentially (exponential phase), reaches a maximum 
size and then declines again. Furthermore, and more 
strikingly, the disadvantage of cooperators can be over- 
come and a transient increase of cooperation can emerge. 
Even though the transient increase is caused by demo- 
graphic fluctuations, it is instructive to examine the de- 
terministic equations first. They not only describe the 
overshoot in the population size well, but also give in- 
sights into the relevant time scales of the dynamics: 



dfX ~ —s{l + px)x{l — x)., (17a) 
dtN ^ (l+px-^\ N (17b) 



The first equation describes the change in the aver- 
age fraction of cooperators. The dynamics occurs on 
the time scale ^ l/s, i.e. the strength of selec- 
tion sets the time-scale of internal evolution. Note that 
df X < always holds and therefore the deterministic ap- 
proximation cannot give rise to any transient increase 
of cooperation. In contrast, the dynamics of the to- 
tal population size is well described deterministically, 
see Fig. [2ja). It resembles the well known equation of 
logistic growth |49| with a frequency-dependent maximal 
population size K{1 + px) (carrying capacity). During 
growth, changes in the population size occur on a time- 
scale Tn ^ 1 + px, c.f. Eq. ( |17b[ ). In the limit of weak 
selection, tn is comparably smaller than the time scale 
Tx, on which selection occurs. This and the frequency 
dependent carrying capacity are the reason for the over- 
shoot: At the beginning the maximal population size is 
given by K{l+pxo). Because cooperators go extinct, the 
size decreases with time. As this reduction is happening 
on a faster time scale than selection, r^r < t^, the popu- 
lation size grows towards a larger carrying capacity, and 
then subsequently drops with decreasing carrying capac- 
ity due to a decline in cooperation. 



B. A Transient Increase of Cooperation 

The stochastic dynamics of the average fraction of co- 
operators (x) is qualitatively different from its determin- 
istic limit. We observe a transient increase in the level of 
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FIG. 2: (Color online) Cooperation in growing populations. 
Temporal development of ensemble averages, (a) The pop- 
ulation size. Starting with A'o = 4, the system grows expo- 
nentially until the carrying capacity is reached. It then falls 
again due to selection and a decreasing carrying capacity, see 
text. The full stochastic solution, gray (red) line, is described 
well by the deterministic approximation, black line, (b) The 
fraction of cooperators. It initially increases due to asymmet- 
ric amplification of fluctuations, and then falls again due to 
selection, see text. The level of cooperation, x, falls below its 
initial value, xq at the cooperation time, tc- The transient 
increase is stronger for larger fluctuations and thus is stronger 
with a smaller initial population size A'o, see gray (colored) 
lines. The deterministic approximation do not account for 
this behavior, cf. black line. Parameters are s = 0.1 and 
p= 10. 



cooperation during a time window {0,tc)- The magni- 
tude of the cooperation time, tc, strongly depends on the 
initial population size A'o, cf. Fig. [2] The origin of this 
transient increase in cooperation is a genuine stochastic 
effect: Demographic fluctuations during the initial phase 
are subsequently asymmetrically amplified by the popu- 
lation dynamics. Heuristically, this can be understood as 
follows; for a detailed mathematical analysis employing 
a van Kampen expansion see the next section. 

For a small initial population size A'o demographic fluc- 
tuations are effectively symmetric, i.e., the occurrence of 
an additional cooperator or an additional free-rider are 
equally likely. However, the consequences of these two 
directions of demographic fluctuations differ strongly: In 
the exponential phase, an additional cooperator ampli- 
fies the growth of the population, while an additional 
free-rider hampers it. In other words, fluctuations to- 
wards more cooperators imply a larger growth rate and 



hence a larger population size. Therefore, those real- 
izations of the stochastic dynamics have a larger weight 
in the ensemble average, Eq. (16b), and enable an in- 



crease in the overall fraction of cooperators. With these 
considerations, a criterion for the transient increase of 
cooperation can be obtained: Demographic fluctuations, 
which are of size ^/N [IT], have to be large enough to 
overcome the selection pressure towards free-riders. This 
can already be inferred from Fig. [2]jb), where curves for 
three different values of the initial population size are 
shown. For the smallest A/q the effect is the strongest 
because fluctuations are large at the beginning. In sum- 
mary, a population bottleneck which corresponds to a 
small initial population size can favor cooperation tran- 
siently. Furthermore, if populations repeatedly undergo 
population bottlenecks, the increase in cooperation can 
be manifested also permanently. 



C. Van Kampen Expansion 

As discussed above the transient increase of coopera- 
tion is caused by fluctuations which are asymmetrically 
amplified. In order to quantify these findings analytically, 
we employ an Omega expansion in the system size ac- 
cording to van Kampen [73] of the master equation ( 15 1. 
For generality, we perform these calculations for arbitrary 
global growth function g{x). The deterministic solutions 
are separated from fluctuations by the following ansatz: 



Nc = nc{t) + Vn^, 



(18) 



c{t) and f{t) correspond to the deterministic solutions, 
as shown below. ^ and fi are fluctuations in the number 
of cooperators and free-riders. The relative strength of 
fluctuations and the deterministic parts are weighted by 
powers of f2 which scales with the current system size. 
For instance, to describe the transient increase which is 
generated by fluctuations at the beginning, is given 
by A'o. Hence, this ansatz accounts for the fact that 

fluctuations scale as \/\/N [73^ Eq. (151 is expanded 

in orders of I/a 



Eg are given by 



Q.. With Eq. (18), the step operators 



E^ 



=1 ± 



1 



1 



E-l± a. 



o(r!3/2). 



(19) 



Employing these and Eqs. ( 18 ) in Eq. ( 15 ) leads to 



+!]"(...) + -!.(...) +0(f^-3/2) 



(20) 
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where terms of the order ^l/K and higher are ne- 
glected. Initially, starting with a small population, 
these higher orders are very small because Vt k, Nq and 
Nq K holds. The orders fi^ and 1/y/Tl depend on 
c, /, s, b, c, g, df^, ^, ^ and are not written 
out in this equation for clarity. By collecting terms 
of order \/Tl and using the identities n = c{t) + f{t) 
and X — c{t) I [c{t) + f (t)] the deterministic equations, 
Eqs. (17 1, are obtained (for K — > oo). Higher orders of 
Eq. (20) lead to a Fokker-Planck equation for 
From this Fokker-Planck equation, differential equations 
for the first and second moments of the fluctuations can 
be obtained. The first moments are given by 



(0 = 



-f- x{l — x)dx- 



x^d. 



9{x)(t>c 



1 



[{l-xf{e)-2x{l-x){S,^^)+x^{^?)] X 



(A) =(1 - xfd.. 

1 



-(0 



=— — x{l — x)dx^^-Y~ 



[{l-xf{e)~2x{l-x){^i,)+x\i?)] X 



(21) 



Note that the second moments only couple at order 
l/^/Tl. Neglecting these higher orders, Eq. (21 1 is lin- 
ear and has an unstable fixed point at (^, n)* = (0, 0). 

Next, we analyze the impact of the second moments 
on the dynamics. Their coupling into Eq. (21) is only 



important for small times, when the first moments are 
still at the initial condition, the unstable fixed point 
(CjM)* — (OjO). Therefore, it is appropriate, to exam- 
ine the second moments for small times, t — > 0. They 
then have the asymptotic form 



dtie 



-Ing^^x, 



(22) 



Due to the inhomogeneity of the differential equations, 
the second moments (^^) and (/^^) immediatly start to 
grow. These non-zero second moments now couple back 
into the first moments, Eqs. ( |21[ ), and push them out 
of the unstable fixed point. To quantify this, the so- 
lution of Eqs. (22) is employed in Eqs. (21). The re- 



sulting equations are solved for small but finite times. 
As the increase of cooperation is caused by fluctuations, 
fluctuations have to establish first. As fixed time we 
here consider the doubling time of the initial popula- 
tion td = l/g{x). Within the time window [0,^^] evo- 
lution is neutral (s <C 5(x))and thus a; = holds. The 
approximation leads to a lower bound for the strength 



of fluctuations. Furthermore, the initial conditions are 
given by (^o) = (mo) = (Co) = (Mo) = 0. If the initially 
generated and asymmetrically enhanced fluctuations are 
large enough to overcome the selection disadvantage, the 
transient increase of cooperation arises. To quantify this, 
the total fraction of cooperators in the system has to be 
examined: 

±u) = (^g) x + i/Vn{i) 

dV""' {Nc+Np) n + l/VnaO + ifJ.)) 

{xn + l/VUm) + (fi) 



For 



dt ^ 



(23) 

> the transient increase of cooperation is 



The condition 



present, j.^^ i.>jiiu.n,±vjii 
s, to the transition line 

dx ln[5(a;)] 



leads, to first order in 



dx9{x) 



n{l/g{x))VL xo n{l/g{x))ng{x) 



(24) 



Here, D, is given by the initial population size Nq. 
For smaller s there is a transient increase in cooper- 
ation, while for larger s the level of cooperation de- 
creases immediately. This resembles the condition for 
neutral evolution, e.g. [151 [52] ; evolution is only neutral 
for sN < const. Thus, only if fluctuations are strong dur- 
ing the initial phase of the dynamics, such that the sys- 
tem behaves neutrally, they are sufficient to overcome the 
selection pressure towards free-riders. The phase bound- 
ary and thereby the strength of the transient increase 
depends on dxg{x)\xQ and g{xQ). Both terms have an- 
tagonistic impacts on the transition line. The reason for 
this behavior is that the initial doubling time, i.e. the 
time during which fluctuations are the most pronounced, 
decreases with increasing g{xQ). The positive enhance- 
ment relies on the growth advantage of more cooperative 
realizations, which depends on dxg{x)\xo at the begin- 
ning. Note, that for non-linear growth functions, where 
dxg{x)\xQ also depends on xq, the transient increase can 
even be reduced by accounting for higher orders. This be- 
havior was also experimentally observed in recent studies 
with microbes, where the growth advantage of coopera- 
tors was tuned [72]. In the next paragraph, we show 
that the calculated phase boundaries match our simula- 
tion results very well for several distinct global growth 
functions. 



D. Phase Diagrams 

In the following we consider how the duration tc of the 
transient increase in cooperation depends on the system 
parameters for the specific global growth function g{x) — 
1 + px, cf. Fig. [3j Then, the transition line between a 
transient increase, tc > 0, and an immediate decrease. 



tc = 0, given by Eq. (24), now reads, 

P 

s — 

nil(l H- pxo) ' 



(25) 
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FIG. 3: (Color online) The transient increase of cooperation 
and its dependence on parameters. Encoded in gray (colored) 
scale, the cooperation time tc is plotted for three different 
pairs of parameters: {A'ojS}, {No,p}, and {xo,s} in (a), (b), 
and (c) respectively. The boundary between the regimes of 
transient increase and immediate decrease are in good agree- 
ment given by Eq. (241, plotted as black lines. In the inset 



of (a), the cooperation time is shown for varying selection 
strength s: tc sharply drops at the boundary. Not varied 
parameters are given by p = 10, Xo — 0.5 in (a); s = 0.05 and 
xo = 0.5 in (b); p = 10, A'o = 6 in (c). 



where nil = 2Ao. For smaller selection strength, s < 
nn(i+pxo) ' asymmetric amplification of fluctuations is 



TiSl(l+p2;o) 



sufficient to overcome the selection disadvantage of coop- 

erators while for larger selection strength, s > — 

free-riders prevail. 

In Fig.jSj we compare this result of the analytical calcu- 
lations with the stochastic simulations. We observe that 
upon increasing the strength of selection, s, which sets 
the advantage of free-riders, the cooperation time tc de- 
creases. In contrast, stronger demographic fluctuations, 
their strength scales as l/v'TVo, prolong the duration of 
the transient increase, i.e., tc increases with decreasing 
A'o, cf. Fig. [sjja). These two antagonistic effects lead to 
a sharp phase boundary between the regimes of transient 
increase {tc > 0) and immediate decrease {tc = 0); see 
inset of Fig. [sjja). Here, the cooperation time steeply 
drops to zero if the strength of selection exceeds a criti- 
cal value. The boundary line is in good agreement with 
Eq. ([25]), cf. black line in Fig.[3];a). 

In Fig. [sjb), the cooperation time is shown for vary- 
ing initial population size A'o and strength of the global 
fitness advantage due to cooperators, p. Now, the phase 



boundary is determined by the interplay between the size 
of demographic fluctuations and its amplification due to 
the global fitness advantage of more cooperative popula- 
tions. A'o has to be small enough for the asymmetric am- 
plification mechanism to be effective. Again, the phase 
boundary is in good agreement with Eq. (24); see solid 
black line in Fig. [Sjjb). 

In Fig. [sj^c), the cooperation time is plotted for vary- 
ing initial cooperator fraction, xq, and selection strength, 
s. We find that the cooperation time decreases with in- 
creasing Xq. Remarkably, for small Xq, the amplification 
mechanism is especially pronounced and therefore able to 
compensate comparably large selection strengths s. This 
is again well described by Eq. (25), see Fig. [3];c) (solid 
black line). The observation is of possible relevance for 
the evolution of cooperation since it allows one a small 
initial fraction of cooperators to proliferate in the popu- 
lation. 

Taken together, our analytical calculations provide a 
mechanistic understanding for the transient increase of 
cooperation and its dependence on the system parame- 
ters s, p, Xq, and A^o- We have quantitatively calculated 
the phase boundary and gained insights into the basic 
nature of the transient increase: First, the probability 
distribution in the cooperator fraction (x) is broadened 
due to neutral evolution; note that Eq. (25 ) resembles the 
condition for neutral evolution [45^ .82^. Second, these 
initially generated fluctuations are asymmetrically am- 
plified and can, therefore, cause an increase in the level 
of cooperation. 



E. The Dormancy Scenario 

Let us now consider the dormancy scenario where the 
ability to reproduce decreases with increasing population 
size. For specificity, we assume the global birth and death 
functions to be given by 



g{x, N) ~ 1 + px — 



N 



and d — 0. 



(26) 



In this scenario individuals do not die but the birth rates 
decrease towards zero as the population size reaches its 
carrying capacity. The relative functions, fs and ws, are 
the same as before; the weakness terms are constant and 



the fitness terms given by Eq. (12). 

To understand the differences in the evolutionary out- 
come, we again study the deterministic rate equations 
first. They are given by 



dtN= [l+px-^]N, 



N 



(27a) 



dtX = -s\^l+px-^—jx{l-x). (27b) 

The equation describing population growth is formally 
identical to the corresponding equation in the balanced 



growth scenario, Eq. (17b). Differences arise because in 
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FIG. 4: (Color online) The dilemma of cooperation in the 
dormancy scenario, (a) The growth dynamics. Initially, the 
small population grows exponentially until growth is stopped, 
cf. light gray (red) line. This behavior is well described by 
the deterministic equation (27a I, see black line. In contrast. 



for the balanced growth scenario, the dynamics continue and, 
due to selection, the population size falls again, see dark gray 
(blue) line, (b) The fraction of cooperators. Equal to the bal- 
anced growth scenario, dark gray (blue) line, there is an ini- 
tial increase of cooperation due to asymmetric amplification 
within the dormancy scenario. Again, this is not described 
by the deterministic approximation, Eq. (27b I. However, in 



contrast to the balanced growth scenario, the higher level of 
selection is latter fixed due to the stop in growth dynamics. 
Parameters are given bys = 0.05, p=10, and A'^o = 4. 



the present case there is mutual feedback between inter- 
nal and population dynamics. This couphng implies that 
both arrest once the population size reaches its carrying 
capacity. In the arrested state there is a relation between 
population size iV* and composition x*: 1+px* = N* /K. 
Thus, the reached stationary state, {x*,N*), depends 
on the initial values xq and A^o- The precise mapping 
depends on the selection strength s. For weak selec- 
tion (small s), the population dynamics is much faster 
than the internal dynamics and hence the population size 
reaches a stationary state while the composition is still 
at its initial value xq, i.e., A^* = K{1 + pxo). In con- 
trast, for strong selection, cooperators go extinct quickly 
with x* = such that the stationary population size 
becomes A^* — K. An example for the deterministic dy- 
namics is shown as a solid black line in Fig. |4] As for 
balanced growth, the deterministic dynamics exhibits a 
strictly monotonous decrease in the cooperator fraction, 
with the difference that now the asymptotic value is ar- 



rested at some finite value. These differences are also 
reflected in the stochastic dynamics, where the asym- 
metric amplification mechanism is acting, cf. Fig. |4j In 
the initial phase of the dynamics, this mechanism affects 
the time evolution of the cooperator fraction in the same 
way as for balanced growth, namely it leads to an initial 
increase of cooperation. Differences in birth and death 



rates, Eq. (13) and Eq. (26), are negligible for small pop- 



ulation size. The arrest of the dynamics only becomes 
effective at later times where an increase in population 
size implies a significantly declining birth rate. As a con- 
sequence even the stochastic dynamics becomes arrested 
such that the initial rise in the cooperator fraction may 
become manifested as a permanent increase. This will 
be the case if the dynamics becomes arrested during the 
time window where the asymmetric amplification mech- 
anism acts; see red line in Fig.Qb). 

In summary, there are now three scenarios for the dy- 
namics, cf. Fig. [5] In addition to the immediate de- 
cline and transient increase there is now also a perma- 
nent increase in the cooperator fraction. The analytical 
expression separating the regimes of transient increase 
and immediate decline still holds, Eq. (25 1, because it is 
due to the same mechanism as before. We did not man- 
age to derive an explicit expression for the transition line 
to permanently increase. However, as the existence of a 
permanent increase in the cooperator fraction depends 
on the asymmetric amplification mechanism, the regime 
of permanent increase is bounded by a hyperbolic line 
beneath the one given by Eq. (24). The latter is a nec- 



essary but not a sufficient condition for the permanent 
increase to occur. 



IV. CONCLUSION 

In this article, we have given a synthesis of evolution- 
ary and population dynamics. This is based on the un- 
derstanding that birth and death events are the driv- 
ing forces underlying changes in the size as well as the 
composition of a population ll(Jl. Both processes are in- 
herently stochastic and inevitably lead to demographic 
fiuctuations whose magnitude depends on the popula- 
tion size. The ensuing stochastic formulation thereby 
naturally accounts for the coupling between internal evo- 
lutionary dynamics and population dynamics. The evo- 
lutionary outcome of the dynamics is determined by the 
interplay between selection pressure and random drift 
caused by demographic fluctuations. Since our approach 
allows to study evolutionary dynamics with varying pop- 
ulation size we can explore ecological situations where the 
relative impact of deterministic and stochastic evolution- 
ary forces change with time. Thereby demographic fluc- 
tuations may lead to a dynamics which is qualitatively 
different from the corresponding deterministic dynamics: 
Beyond creating a broad distribution in size and compo- 
sition, the coupling can strongly distort the distribution 
and thus strongly influence average values. For the public 
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5 10 15 20 25 30 35 

(initial group size) Nq 

FIG. 5: (Color online) The transient increase of coopera- 
tion for the dormancy scenario. The cooperation time tc de- 
pending on the initial population size, A'^o, and the strength 
of selection, s. The condition for a transient increase of co- 
operation to occur is still given by Eq. (241, black line. In 



addition, due to the stop in growth dynamics, there is an ad- 
ditional regime, where the increase becomes permanent, dark 
gray (dark blue) area. The permanent increase is also shown 
in the inset, where the cooperation time is shown for varying 
strength of selection. If, for a given initial population size, se- 
lection is sufficiently slow compared to fixation of the growth 
dynamics, the increase of cooperation becomes permanent. 
Parameters are given by p = 10, xq — 0.5. 



good scenario, discussed in this paper, this corresponds 
to an asymmetric amplification mechanism which yields 
a transient increase in the level of cooperation. 

In the absence of a coupling between internal evolution 
and population dynamics, the impact of population size 
on the internal evolutionary dynamics reduces to a mod- 
ulation in the strength of demographic fluctuations. If, in 
addition, the deterministic population dynamics exhibits 
a strongly attractive fixed point at a finite population 
size, our model maps to a standard description of evolu- 
tionary dynamics, i.e. the Moran process. 

The general observations made for the coupled stochas- 
tic dynamics are exemplified by the dilemma of coopera- 
tion in growing populations. Here, fluctuations in combi- 
nation with growth lead to a transient increase of cooper- 



ation. The origin of this increase is the asymmetric am- 
plification of fiuctuations. As the presence of cooperators 
increases the growth rates, fluctuations towards those are 
enhanced. Therefore growth dynamics cannot be ignored 
but can be an essential part in evolution. Furthermore, 
the details of the growth dynamics can be crucial in de- 
termining the evolutionary outcome. As we have consid- 
ered for the dilemma of cooperation and two extremes 
of microbial growth dynamics, cooperation can either in- 
crease only transiently or the higher level can even fixate 
due to dormancy. Our analytical derived transition line 
provides the same sufficient condition for the transient 
increase in both scenarios. Furthermore, the same line is 
also a necessary condition for the permanent increase for 
the dormancy scenario. In actual populations, both sce- 
narios are present with a fraction of 20% - 80% dormant 
bacteria [93]. While the transient increase does not de- 
pend on this fraction, the permanent increase is smaller 
than for purely dormant bacteria. The discussed scenar- 
ios for the increase of cooperation, rely on demographic 
fiuctuations which are especially pronounced during pop- 
ulation bottlenecks. Such bottlenecks may be caused by 
seasonal changes of the environment, migration into new 
habitats and range expansion, e.g. [iSl l^5H5^ . In addi- 
tion, if the permanent increase is not present, repeated 
bottlenecks provoking regular occurring growth phases 
can favor cooperative behavior by stabilizing a former 
transient increase. This becomes especially important in 
the context of biofilms where population structure and in- 
volved restructuring mechanisms can drastically change 
evolutionary outcome [SSI IZIl El] ■ 
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